A spectral approach to the tunneling decay of two interacting bosons 



Stefan Hunn/ Klaus Zimniermann,^ Moritz Hiller/'^ and Andreas Buchleitner^ 

^ Physikalisches Institut, Albert-Ludwigs-Universitat Freiburg, 
Hermann-Herder-Str. 3, 79104 Freiburg, Germany 
^Institute for Theoretical Physics, Vienna University of Technology, 
Wiedner Hauptstrafle 8-10/136, IO4O Vienna, Austria 
(Dated: March 5, 2013) 

We study the full-fledged microscopic dynamics of two interacting, ultracold bosons in a one- 
dimensional double-well potential, through the numerically exact diagonalization of the many-body 
Hamiltonian. With the particles initially prepared in the left well, we increase the width of the 
right well in subsequent trap realizations and witness how the tunneling oscillations evolve into 
particle loss. In this closed system, we analyze the spectral signatures of single- and two-particle 
tunneling for the entire range of repulsive interactions. We conclude that for comparable widths 
of the two wells, pair-wise tunneling of the bosons may be realized for specific system parameters. 
In contrast, the decay process (corresponding to a broad right well) is dominated by uncorrelated 
single-particle decay. 



PACS numbers: 03.75.Lm, 03.65.Xp, 05.30.Jp 

I. INTRODUCTION 

One of the hallmarks of quantum mechanics is tun- 
neling through energy barriers that classically cannot be 
overcome. It is at the heart of fundamental dynamical 
phenomena like Josephson oscillations [I] , the spatial flip- 
ping of the nitrogen atom in the ammonium molecule [2] 
or quantum lattice dynamics [3] . Tunneling is also one of 
the major mechanisms that provoke the loss of particles 
from quantum systems [4]. The latter occurs whenever 
a quantum system is coupled via a potential barrier to 
a continuum, i.e., to a set of asymptotically free states. 
Such particle loss is omnipresent in nature, e.g., in ra- 
dioactive decay or in the autoionization of excited atomic 
states. 

While the tunneling and loss dynamics are thoroughly 
understood for single particles [5., much less is known 
about the generic case of interacting many-body systems. 
Furthermore, in the vast majority of physical realizations, 
the inter-particle interactions are immanent and, if at all, 
hard to manipulate. An archetypical example is the inter- 
electronic Coulomb interaction which plays a crucial role 
in the non-sequential double ionization of helium [6], or 
in the sequential double ionization of argon [7]. 

In this respect, the advent of ultracold atoms in opti- 
cal potentials prepared the stage for a new generation of 
experiments in which all relevant control parameters can 
be controlled essentially at will, ranging from the poten- 
tial landscape [5HTU]. over the initial state (and number) 
of atoms [TT] , to the inter-particle interaction, which can 
be tuned over several orders of magnitude from attrac- 
tive to repulsive [12l [13] . Consequently, several authors 
have studied the decay of particles from a Bose-Einstein 
condensate (for recent experiments, see, [TH [TS]). The 
two major theoretical tools used in this context are the 
mean- field (Gross-Pitaevskii) description [T6lf22] . and ef- 
fective Hamiltonian or master equation approaches |231 - 



[27], to name a few of them. While both methods can 
yield valuable insight, they face fundamental limitations: 
The mean-field treatment (applicable in the limit of large 
particle numbers) is based on the assumption that the 
condensate is in a coherent state and breaks down p8]- 
I31j when the inter-particle interaction leads to consid- 
erable dephasing 133] . The master equation, on the 
other hand, is typically introduced in an ad hoc procedure 
without a rigorous justification and few authors went be- 
yond the standard Markovian treatment |31H3B]. 

In contrast, the microscopic decay process of interact- 
ing particles only begins to be explored on a fundamen- 
tal level. Recent experiments reported the interaction- 
dependent loss of the first of a few atoms from a one- 
dimensional trap |37j . and the corresponding escape 
rate was successfully described by a quasi-particle wave- 
function approach [BS]. The theoretical works that de- 
scribe the full-fledged decay dynamics rely on the propa- 
gation of the time-dependent Schrodinger equation either 
by direct integration (3^ , by combining a matrix-product 
state approach with Bose-Hubbard-like chains [30] : or via 
multi-conflgurational Hartree-Fock methods [U [32] • 

In our present work, we take a complementary ap- 
proach to tunneling decay and study the numerically ex- 
act many-body dynamics of two interacting bosons which 
are initially prepared in the left site of a double- well po- 
tential. In subsequent realizations of the trapping po- 
tential, we gradually increase the extension of the right 
well which — in the limit of large widths — mimics the 
unconfined configuration space (to which the particles es- 
cape), by a dense quasi-continuum of states. Of course, 
the finite extension of the broad well defines a maximal 
observation time before the particles hit the rightmost 
boundary. This limitation is, however, more than com- 
pensated by two major benefits: On the one hand, we can 
directly monitor the transition from the regime of tun- 
neling oscillations between the two wells to the regime of 
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tunneling decay from the left well, in a single setup. On 
the other hand, we gain complete access to all spectral 
quantities of the system (e.g. its eigenvalues, eigenfunc- 
tions, and the density of states) thus providing insight to 
the inner workings of many-body tunneling. 

The key problem is to expose the physical nature of the 
many-body decay process: Given, for example, the quan- 
tum mechanical correlations naturally present in the ini- 
tial state of two interacting bosons, do the particles leave 
the left well as a pair or do they tunnel independently? 
That is, depending on the strength of the interatomic in- 
teractions, will we observe predominantly correlated or 
uncorrelated processes? We would as well like to know 
which quasi-continuum states actually support the de- 
cay process, and how interaction energy is converted into 
kinetic energy, during the tunneling decay. To answer 
these questions, we shall investigate two-body as well as 
reduced single-particle quantities. 

The narrative of this work is as follows: In the next 
section we introduce our model, the double-well poten- 
tial, and those quantities used throughout the paper to 
study particle loss. Sec. Ill is devoted to the tunneling 



decay of a single particle and sets the frame of reference 
for the two-body problem. From a spectral perspective, 
we analyze the transition from tunneling oscillations to 
tunneling decay as the width of the right well is increased, 
and will see that the participation number of the initial 
state in the eigenstates of the double well is an excellent 
figure of merit to characterize this transition. We briefly 
revise single-particle decay theory, and elaborate on the 
validity of our quasi-continuum approach. 

In Sec. |IV| we introduce the initial state of the dynami- 
cal two-boson evolution, and define the principal physical 
observables. Our central results are presented starting 
from Sec. [V] where we analyze the spectral properties of 
the tunneling process of two interacting bosons as the 
width of the right well is gradually increased. Specif- 
ically, we identify those particle configurations in the 
quasi-continuum states which support single- and two- 
particle tunneling, respectively, and calculate the associ- 
ated density of states. This procedure parallels the anal- 
ysis of Sec. |III[ and the differences to the single-particle 



case are exposed. In Sec. VI the predictions of our spec- 
tral analysis are compared to the numerically obtained 
time evolution. We develop a complete picture of the 
decay process by means of two- and single-particle quan- 
tities and compare our results to related studies on the 
time-resolved decay dynamics [39l [HI [42] . We conclude 
with a summary and discussion in Sec. |VII[ 



II. DOUBLE- WELL POTENTIAL AS A MODEL 
FOR TUNNELING DECAY 

We now prepare our tools to describe the tunneling 
dynamics of two ultracold, interacting bosons in a one- 
dimensional double-well potential. For x and y the re- 
spective coordinates of the (indistinguishable and spin- 



a) 



V{x) 



® ® 



Vn 





/ 


) 






y 












y 





1 — ► 



r 



X 



b) 



y 



(2) 






(2) 



FIG. 1. (Color online) Top: (a) Schematic sketch of the 
single-particle potential V{x) as defined in Eq. Q. The left 
well of width £ (interval I) is coupled via a potential barrier 
of width b and height Vb (interval II) to the right well of 
width r (interval III). Bottom: (b) Potential landscape of (ij) 
in the two-boson x — y configuration space. Two particles 
confined in the left well correspond to region (1) (black), one 
particle in either well corresponds to the two regions (2) (red 
or light gray), and both particles in the right well corresponds 
to region (3) (blue or dark gray). The tunneling barrier is 
depicted in white. The dashed white line denotes the diagonal 
{x = y) at which the contact interaction U between the two 
particles arises. 



less) bosons, the two-particle Hamiltonian reads 

Htpix, y) = Hsp{x) + Hsp{y) + W{x, y) , (1) 

where Hgp is the single-particle Hamiltonian and W{x^y) 
accounts for the inter-particle interaction (see also Ap- 
pendix The former is given by 



Hsp{x) = - 



da;2 



V{x) 



(2) 



where we have fixed the particle mass m = 1/2 and h = 
1. The corresponding single-particle potential Vix) is 
depicted in Fig. [l^): the left well of width I (interval I) 
is coupled via a potential barrier of width 6 and height 
Vb (interval II) to a second well of width r (interval III) , 
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The interaction term W{x, y) is assumed to result pre- 
dominantly from s-wave scattering [43], due to the ex- 
tremely low temperature of the bosons. This leads to a 
(5-like interaction between the particles, i.e., W{x,y) — 
U5{x — y) with U — 2hujradas Here, as is the s- 

wave scattering length, which can be tuned experimen- 
tally over several orders of magnitude, from attractive to 
repulsive [12^ il3j , and ujrad denotes the frequency of the 
harmonic radial trapping potential which acts transver- 
sally to the direction of For a cigar-shaped trap (i.e., 
a small longitudinal trap frequency uJiong ^ ^rad) we can 
effectively treat the potential as being one dimensional, 
since — due to the low total particle energy — only the 
transversal ground state is occupied. To define the units 
(of energy etc.), we introduce a reference length scale 
C With the above convention on h and m, energy is thus 
measured in units of time in units of and length 

in units of C. The interaction strength U carries the di- 
mension energy times length is thus measured in units of 
l/£. 

To investigate tunneling decay (or, synonymously, par- 
ticle loss) we consider the following dynamical scenario: 
Initially, the particles are prepared in the many-body 
ground state of the isolated left well, corresponding to an 
infinitely high tunneling barrier Vq = oo (or b — r — 0). 
Then, the barrier height is instantaneously reduced to a 
finite value Vq — 0.1. We study the ensuing time evo- 
lution for different trap configurations, at fixed widths 
i = 51 and 6 = 2, but for ever broader right wells. The 
values for b and Vq are chosen such that significant tun- 
neling takes place within comparably short times. At the 
same time, this choice ensures that the energy of the ini- 
tial state is always considerably smaller than the height 
Vq of the barrier, such that quantum tunneling is the only 
way to leave the left well. 

In this way, we define a closed setup, in which the left 
well serves as the system while the broad right well pos- 
sesses a large density of states and can be regarded as the 
environment to which the bosons escape. Throughout 
this work, environment is to be understood as a quasi- 
continuum of (two-boson) states, i.e., as the unconfined 
configuration space. The rational behind this approach 
is twofold: On the one hand, we can follow the crossover 
from tunneling oscillations (expected for r w £) to parti- 
cle loss — realized in the limit of a large width r ^ £. On 
the other hand, the corresponding Hamiltonians ([T]), (2| 
can be numerically exactly diagonalized (see Appendix 
|A] for details). The price we pay is a finite observa- 
tion time, defined by the condition that no reflections 
from the hard wall boundary of the right well must oc- 
cur. We are, however, rewarded with the direct access 
to all spectral quantities of the problem. This represents 
a complementary approach to the direct solution of the 
time-dependent Schrodinger equation , and yields 

fundamental insight as to which quasi-continuum states 
actually support the decay. This information should also 
advance the derivation of a microscopic master equation 
and we shall shortly return to this at the very end of this 



article. 

The major theme of our study is under which circum- 
stances and to which extent the bosons tunnel as indi- 
vidual particles or in a correlated manner, i.e., as a pair. 
To facilitate the forthcoming discussion, we visualize the 
potential landscape of ([T]) in the x-y-plane, see Fig.[l|D). 
The white regions denote the tunneling barrier of height 
Vb, and the dashed line represents the diagonal {x — y) 
at which the contact interaction U arises. Region (1) 
corresponds to both particles being in the left well, while 
regions (2) indicate that one boson is in either well. Fi- 
nally, region (3) corresponds to both particles being lo- 
cated in the right well. In this representation, uncorre- 
lated tunneling of the two bosons would manifest in a 
transition (1) — > (2) — ^ (3) while correlated tunneling (of 
a pair of bosons) would correspond to a direct transition 
(1) (3). Further down, in Eqs.( [23p5| , we shall define 
the associated probabilities with which the two-particle 
wave function populates the corresponding regions. 

As an aside, we point out that the above considered 
transitions can be regarded as sequential [(1) (2) 
(3)] or non-sequential [(1) — > (3)] double ionization of 
many-electron atoms. The continuing interest in the He- 
lium problem is driven by the same question as in the 
present study: That is, to identify the role of the inter- 
electronic interactions and correlations in the ionization 
process (for recent publications, see, e.g. 



III. SINGLE-PARTICLE CASE 

As a preliminary step toward the tunneling of two in- 
teracting bosons, we revise the single-particle scenario of 
Hamiltonian ^ and investigate the spectral signatures 
of the transition from tunneling (r £) to particle loss 
(r ^ £). This allows us to introduce — in the familiar 
single-particle context — one of our main spectral tools, 
the participation ratio (PR) of the initial state in the 
Hamiltonian's eigenbasis, i.e., the number of eigenstates 
that mediate the dynamics. The analysis of the PR con- 
stitutes the first part of the present section and will prove 
beneficial once we turn to the two-particle dynamics. At 
the end of this section, we demonstrate the validity of 
the quasi-continuum approach, by comparing our results 
of the single-particle spectrum to the time evolution of 
the initial state and to analytical predictions. 



A. Initial state 

Throughout this contribution, we focus on the ground 
state of the isolated left well as the initial state \tp{t = 0)) 
of the dynamics. For a single particle in a box potential 
of width £, the configuration space representation of the 
eigenstates \Xn) are sinusoidal [15], i.e.. 
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where kn^"^ = irn/i denotes the single-particle momen- 
tum and the corresponding eigenenergies are 



(5) 



Thus, with the numbers from above the ground-state en- 
ergy e^^^-* = 3.79 X 10^"^ is well below the barrier height 
Vq = 0.1. Here and in the following, the super-script (sp) 
indicates single-particle quantities. If not stated other- 
wise, the quantities kn^^ and e^*^'' pertain to the left well, 
hence we drop the index £. 

Since the |xl) are typically not eigenstates of the 
double-well potential (also referred to as the total sys- 
tem in the following), we generally observe non-trivial 
dynamics once we set \ip{0)) = \Xn)- Since numerically 
exact diagonalization is at the heart of our approach, we 
obtain the time-evolved state by spectral decomposition. 
That is, we expand \ip{0)) in terms of the eigenstates 
{\Ek'P^)} of the total system 



= (i?(^P)|^(0)), 



and use the familiar expression 
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FIG. 2. (Color online) Single-particle case: The participation 
ratio PR (magenta), and the weighted participation ratios 
PKe (black) and PR,, (blue) for different widths r of the right 
well [see Eqs. (|8|, ([lO|]. The dashed green line denotes a linear 
fit with slope a — 1.02x 10"'^. The participation ratios exhibit 
resonances which increasingly overlap as r — >■ oo. Inset: The 
absolute value squared (black circles) |c„p of the expansion 
coefficients ([6| of the initial state, versus the level index n, 
for a broad right well of width r = 3000. The blue solid line 
represents a Lorentzian fit of width F = 1.3. 



B. Participation ratio I: From tunneling 
oscillations to tunneling decay 

A robust and intuitive spectral measure to character- 
ize the quantum dynamics is the participation ratio (or 
participation number) 



mm 



El 



(8) 



which represents the number of eigenstates lEn^"^) that 
significantly contribute to the initial state \ip{0)). The 
participation ratio varies from unity -when |V'(0)) coin- 
cides with one eigenstate- to N, for \ip{0)) an equally- 
weighted superposition of all eigenstates. Here, N is the 
dimension of the associated Hilbert space which is always 
finite in our numerical study (see Appendix [A|). 

Yet, as defined in Eq. ([s]), the participation ratio does 
not take into account the spatial density distribution of 
the \En^''), i.e., the PR cannot resolve which part of con- 
figuration space is populated by the contributing eigen- 
states. Such a distinction, however, is highly desirable 
since it can provide valuable insight into the systems dy- 
namics, as we shall see below. To achieve this, we employ 



the probability Pg to find the particle described by \ip{t)) 
m the left welH 
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With the help of ^ we can define weighted versions of 
the PR as follows: 



(sp)\\2 



EJc„pp.,,(|i?^)))]2j 

X V|c„pP,,,(|£^))). 



(10) 



The numerator in the first term of ( 10 ) weighs each coef- 
ficient I c„ p with Pi{\E„)) [Pri\E„) ) ]7i • e . , with the prob- 
ability that a particle in the corresponding eigenstate 
\En) is found in the left (right) well, while the denom- 
inator fixes the range of the first term in ( 10 ) between 
one and N. The last term in Eq. ( 10 ) represents an 



additional weighting which we explain in the discussion 
around Footnote [i] Put differently, PR/(|?A)) [PR^d?/'))] 
provides a measure for the number of those eigenstates 
that participate in the time evolution and have a fraction 
in the left [right] region of configuration space. 



In the two-boson dynamics studied from Sec. |IV| onwards, we ^ Note that the probability density within the barrier is negligible, 
expand the (two-body) initial state l |17| in terms of the two- Hence, Pr can be regarded as the probability to find the particle 

particle eigenstates of the total system. in the right well. 
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With these quantities at hand, we now discuss the tun- 
nehng dynamics of a particle initially prepared in the 
ground state \xi) of the isolated left well, in Fig. [2j for 
variable width r of the right well. 

Consider first the familiar case of a symmetric double- 
well potential [r — I = 51). The PR takes the value two 
which can be easily understood by the following argu- 
ment: In the symmetric case, the ground states |xi)i IXi) 
of the isolated left and right well are in resonance. Thus, 
the eigenstates \En^^) of the total system are the sym- 
metric and antisymmetric superpositions of the isolated 
wells' eigenfunctions, such that the two lowest states read 
~ (Ixi) ± IxI))/a/2- Accordingly, the initial state 
IXi) is a superposition of these two eigenstates and hence 
the PR equals twoj^ As a result, the associated tunneling 
dynamics displays perfect Rabi oscillations with a period 
inversely proportional to the energy splitting E^^'^ — E'j'*^-' 
(not shown). How does the weighted participation ratio 
PRf (black curve) behave? Said the above it is clear 
that two eigenstates contribute to |'!/'(0)) which (par- 
tially) populate the left well, thus the first term in ( 10 ) 
equals two. Since each of these states amount to half the 
norm of the initial state and Pi{\Ei)) = Pi{\E2)) = 1/2, 
the second term equals 1/2. In total we have PR^ = 1 
and, due to parity symmetry, also PR^ = 1 (blue curve). 

As r increases to, say, 75, the PR approaches one, in- 
dicating that Ixi) is approximately an eigenstate of the 
double- well system. The physical origin is that at r = 75, 
the eigenstates of the isolated right well are maximally 
detuned from resonance with |xi)- Accordingly, PR,, 
tends to zercj^ and PR^ « 1. The associated dynam- 
ics is trivial: the particle is persistently trapped in the 
left well and only small amounts of probability leak out. 

Once r is further increased, the PR exhibits a sequence 
of equally spaced peaks. These occur whenever an eigen- 
state of the isolated right well |x5i) comes into resonance 
with the ground state of the isolated left well |xi), i-e., 
for r sa ni, n — 1,2,... . Accordingly, we will refer to 
the maxima as resonances. At these values of r, one ob- 
serves Rabi oscillations while, in between two peaks, the 
tunneling is suppressed}^ 

The situation drastically changes once the width r is 
substantially increased. Then, the oscillations in the PR 
are reduced, the PR and the PR^ become larger than one. 



^ Throughout the paper we consider the weak couphng regime (see 
Appendix [d||, otherwise more states could contribute to \E^^2^), 
even in the case of the symmetric double well. 

^ In this case all coefficients {|cnpfr{|£'n))} are approximately 
zero. However, due to the normalization via the denominator 
in Eq.| |10| l, the PR^ could become very large. This artificial 
increase is prevented by the additional weighting with the last 
term in Eq.|[Io|. 

^ We note that a similar situation arises in lattices systems sub- 
ject to an additional tilt: At appropriately tuned values of the 
lattice tilt, neighboring sites become resonant, what leads to the 
resonantly-enhanced tunneling (see, e.g. |49p. 



for all r, while PRf decreases. In this regime, the widths 
of the resonances (between eigenstates of the left and 
the right well) are larger than their respective spacing, 
which is simply given by the difference of two consecutive 
eigenenergies tn^'' of the right wellj^ In our setup, this 
transition takes place around r « 1000, that is, r > 1000 
denotes the regime of overlapping resonances, where the 
density of states in the right well is large enough such 
that many states contribute to the initial state, and a 
quasi- continuum is formed; a prerequisite to observe loss 
dynamics rather than Rabi oscillations. Regarding the 
weighted participation ratios in the limit of large r, we 
observe that PR^ saturates at a finite value|^ At the same 
time, PRr approaches PR, since the eigenstates I-EI"^') 
that contribute to the initial state |'0(O)) have an ever 
increasing overlap with the right well. 



C. Participation ratio II: Tunneling decay rate 

Before we turn to the corresponding time evolution, 
let us formulate the expected particle loss quantitatively. 
According to Wigner's theory of decaying systems (see, 
e.g. [5^), the coefficients |c„p should follow a Lorentzian 
distribution of width 7 in energy space [see also the 
schematic illustration in Fig. [7^)] leading to an exponen- 
tial decay behavior in the time domain with rate 7]^ For 
large enough r and sufficiently weak coupling between the 
two wells (as in our case), the mean level spacing of the 
total system A„ = E^+i — En^^ can be safely assumed as 
constant within the extension of the Lorentzian. In this 
way, we can consider the |c„p as a function of the level 

index n rather than of the energy En^\ what is much 
closer to the definition of the PR. In the inset of Fig. [2] 
we plot the |c„p versus n, together with a Lorentzian fit 
(solid line) of width F = 1.3 for r = 3000. The agreement 
with the theoretical prediction is almost perfect and it re- 
mains to connect F to the decay rate 7 and, ultimately, 
to the participation ratio PR. The first step is a simple 



® We note that the width of the resonances increases monotonically 
with r. 

We note in passing that the asymptotic value of PR^ can be 
derived by assuming that, within the energy window in which 
the coefficients |c„p are large, the PR^dE^}) follow the same 
Lorentzian distribution as the coefficients |c„p (see following 
subsection). We obtain an asymptotic value of 16/40 = 0.4 which 
agrees well with our numerics. We numerically confirmed that 
this assumption on the distribution of the PR^ (lEn)) is indeed 
justified for r — 00 (see Appendix 
* We note that this only holds for the initial state \il>(0)) = \Xn) 
an eigenstate of the isolated left well (see AppendixjDji. Further- 
more, we stress that the exponential decay is realized on inter- 
mediate time scales, while deviations thereof naturally occur for 
ultra-short and ultra-long times and have even been experimen- 
tally observed [51] [52]. 



6 



0.1 





' 1 


numerics: 

o o fit 

o o analytics 


- \ IX2> 


IXi> 




, <^ 1 


1 


1 



5000 



t 



10000 



15000 



FIG. 3. (Color online) Single-particle decay: Semi- 

logarithmic plot of the probability Pt{t) to find the boson 
in the left well versus time t, Eq. The black (red) line 
corresponds to an initial condition defined by the ground state 
IXi) (first excited state [xi)) of the isolated left well (jij). Cir- 
cles denote exponential fits [Eqs. (15l,(16l] while diamonds 



represent the analytical result (see text and Appendix \B\ 



proportionality relation: 



where 



r 1 
2-K 4e 



(11) 



(12) 



is the single-particle density of states in the quasi- 
continuum, evaluated at the energy E = t^^^ of the ini- 
tial state. In a second step, we calculate the PR for per- 
fectly Lorentzian distributed expansion coefRcients |c„p 
to be PR(|'0(O))) — ttF. Hence, the loss rate becomes 
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PR(|^(0))) ^ .fjly) 
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It is physical to assume that for large enough r, the rate 
7 of the decay should not depend anymore on the width 
r of the right well. We thus conclude from ( [l3| that the 
participation ratio grows linearly with r, i.e., PR/r — >■ 
a = const, for r — > oo. The expected linear increase is 
confirmed by Figj2]where the dashed green line represents 
a linear fit on the PR. From the slope a = 1.02 x 10"'^ 
we obtain the loss rate 



located in the left well [see Eq. ([9])], as the relevant ob- 
servable. The time evolution is obtained from Eq. JTI) 
by diagonalization of the single-particle Hamiltonian 
In Fig. [sj we plot Pe{t) versus t, for \tp{0)) — |xT/- 
We observe an exponential decay, as expected from the 
Lorentzian behavior in the energy domain. The dynam- 
ical loss rate is obtained by an exponential fit to be 



ry['P^ = 1.16 X 10' 



(15) 



and thus found in good agreement with the spectrally 
determined rate ( 14 ) . In anticipation of the two-particle 



dynamics discussed further down, we also consider the 
decay of the first excited state 1X2)1 which will play a 
role when the inter-particle interaction leads to fermion- 
ization of the bosons. The decay of \ip{t = 0)) = Ixi) 
is plotted as a red line in the same figure, with a fitted 
decay rate 



7^'^^ = 8.5 X 10-4 . 



(16) 



Due to its larger energy, Ixi) effectively feels a lower bar- 
rier than Ixi); hence 72 '^^^ is considerably larger as com- 
pared to the rate 7^'''''. 

The question naturally arises whether the results de- 
rived with the above quasi-continuum approach agree 
with the actual loss physics to a perfect continuum, as- 
sociated with a right well that stretches out to infinity. 
Therefore, we analytically solve the related scattering 
problem as a benchmark. We obtain complex eigenen- 
ergies with imaginary parts which be interpreted as loss 
rates (cf. Appendix|B|. The induced decay is represented 
in Fig. |3] by diamonds, and agrees well with the quasi- 
continuum results. We thus conclude that our quasi- 
continuum approach is indeed the appropriate tool for 
our study of the temporal behavior on not too long time 
scales. 

In the familiar single-particle case, we thus successfully 
used the participation ratio to observe the transition from 
tunneling oscillations (isolated resonances) to tunneling 
decay (strongly overlapping resonances) . In the next sec- 
tions, we will generalize this concept to the case of two 
interacting bosons. 



IV. TWO-PARTICLE CASE: INITIAL STATE 
AND PHYSICAL OBSERVABLES 



Ground State of two bosons in a box potential 
as initial state 



e^^'P^ = 1.25 X 10^4 



(14) 



D. Time evolution and analytical treatment 

We now turn to the particle loss dynamics, with the 
probability Pe{t) = P({\ip{t))) that the particle is still 



We now turn to the case of two interacting particles. 
As a first step, we define the initial state of the system, 
i.e., an eigenstate of two bosons in the isolated left well 
of width ^ = 51. The latter can be obtained analytically 
by means of a Bethe ansatz |53| from the corresponding 
Hamiltonian ([T]) with r = = b. More precisely, the n-th 
eigenfunction V'fci „fc2 „ (2;, y) is completely characterized 
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FIG. 4. Various observables of two interacting bosons in a single-well potential of width I = 51: (a) The ground state momenta 
(black) and fc2,i (red) [see (|18[)], in units of the inverse width 1/^ of the well, versus the inter-particle interaction strength 
U. The dashed horizontal lines denote the asymptotic values k£ = vr, and 2n, respectively, (b) Lowest two-particle energy 
levels versus U, obtained fro m th e exact numerical diagonalization of Eq. ([T]). (c) Diagonal part pred{x,x) of the reduced 
single-particle density matrix (201 of the respective two-particle ground state, versus position x, for various values of U. The 
repulsive interaction reduces the probability to find a boson at the center of the trap, at x = 1/2. 



by two wave vectors fci^„ and fc2,n and reads, for y < x [541 155) : 



'0fci,„fc2,„ (2;, y) oc [Ai exp (iki^nx) - A2 exp(-ifci^„a;)] sin(fc2,ny) ^^^^ 
+ [A3 exp {ik2^nx) - exp(-iA:2^„a;)] sin(A:i^„y) . 



The amplitudes Ai{ki „, k2.n) are functions of the two- inter-particle interaction strength U and represent the 
particle wave vectors ki^n which, in turn, depend on the solutions of the coupled equations: 



ki^ni = "n-iTT + arctan 
k2.ni — n27T — arctan 



U 



kl.n ~ k2.n 
U 

kl n ^2,71 



arctan 



arctan 



U 



U 

kl,n + fc2,n 



(18) 



For /7 = 0, the equations decouple and we immediately 
recover the single-particle result in which the two positive 
integers ni^2 are simply the quantum numbers of a single 
particle in a box of size £. In the general case [/ ^ 0, the 
ki^n can be regarded as the momenta of two bosons with 
corresponding two-particle energy 



(19) 



As in the single-particle case, the quantities ki^n and e„ 
pertain to the left well (if not stated otherwise). 

Since we are dealing with ultracold atoms, we are 
predominately interested in the initial state being the 
ground state [i.e., ni = ^2 = 1 in Eq. (18)]. In the re- 
mainder of this subsection, we recall its properties as a 
function of the interaction strength U. We character- 
ize the ground state by the momenta fc^^i and the diago- 
nal part of the associated single-particle reduced density- 



matrix 

Pred{x,x',t) ^ I dy '4!*{x,y,t)'ip{x' ,y,t) . (20) 



The latter quantity determines all single-particle quanti- 
ties as, for example, the number of bosons in a specific 
region of configuration space [compare Eq. (26)]. 

In Fig. |4^), we plot the momenta fci^i and A;2,i of the 
ground state, in units of the inverse width l/£, versus 
the inter-particle interaction strength U. For U = 0, the 
bosons are independent of each other, and both momenta 
take the single-particle value ki^ii = ir. Accordingly, the 
two-particle ground-state wave function at J7 = is a 



® In this section, ip{x,y,t') = ipf^-^ ^^.j ^{x,y) while we consider 
time-dependent wave functions further down. 
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product of single-particle wave functions, 

'^kiTklii^^y) = ^i(^) sin(fcj'*^^x) sin(fc^''''^y) , 

(21) 

with a two-particle ground-state energy ti{U = 0) = 

26^**^^ given by twice the single-particle ground state en- 
ergy (5]). Since both bosons occupy the same single- 
particle state, also the diagonal of the corresponding re- 
duced single-particle density matrix [panel c)] is perfectly 
sinusoidal. 

As U increases, the particles repel each other and the 
two-particle energy grows, until it saturates at the energy 
of the first excited state ^{x,y): this is evident 

from panel (b) where we plot the parametric evolution 
of the lowest-lying energy levels e„, obtained from the 
exact diagonalization of Eq. ([T]), with r — = b and 
I = 5l{^ In the reduced single-particle density matrix, 
this repulsion between the bosons is manifest in a reduced 
probability to detect a boson at the center of the box, i.e. 
a.t X — £/2, see panel c). 

This observation can immediately be understood in 
terms of the wave vectors: one finds that the (larger) mo- 
mentum ki I grows with U, and asymptotically reaches 
ki^il = 27r, while fc2,i first slightly decreases and asymp- 
totically approaches ki^il = n The indistinguisha- 
bility of the particles aside, one may imagine that for 
U ^ oo one boson is in the single-particle ground state 
IXi), while the other one is in the single-particle first 
excited state |x2)i ^ situation also referred to as the 
fermionization or Tonks-Girardeau limit [56^ of a repul- 
sively interacting Bose gas. In this limit, the wave func- 
tion exactly vanishes on the diagonal where the (contact) 
interaction arises and is given by 

^k^iTk^li^'V) °^ 4s [sin(fc^"''^y) sm{k["P^x)- 

v2 L (22) 

sin(A;2''^''a;) sin{k['^^^y) , 



for y < X, while its energy ^^{U — od) — e^''*''' + e 
is the sum of the single-particle ground and first-excited 
state energies. Of course, V'fei ifc2.i (^^j 2/) always retains 
its bosonic symmetry with respect to exchange of x and 
y. Thus, the ground state wave function of the (interact- 
ing) bosonic and of the corresponding (interaction-free) 
fermionic system coincide in modulus, but have opposite 
sign for X > y. 



(sp) 
2 



B. Physical observables 

In our analysis of tunneling decay, we initially prepare 
the ground state \^{t = 0)) = I'f/'fci ifea.i)' ^'i- ^^'^^^ of 



two bosons in the isolated left well (Vq = oo) at a given 
interaction strength ?7, and monitor the ensuing dynam- 
ics after the potential barrier is instantaneously set to 
Vo = 0.1. To distinguish uncorrelated single-particle tun- 
neling from correlated two-particle tunneling, we define 
-according to Fig. fib)- the following observables: The 
probability Pi{\^p{tJ)) that both particles are in the left 
well [i.e. in region (1)], the probability P2(|V'(i))) that 
one boson is in either one of the wells [i.e. in region (2)], 
and the probability P3{\ip{t))) that both bosons are 
located in the right well [i.e. in region (3)], 

Piiim))^ J^dx jdy\i;{x,y,t)\\ (23) 

P2(|V'(i))) = 2 /dx / dy\^{x,y,t)W (24) 
Ji Jiii 



dx 



III 



dy\4,{x,y,t)\\ (25) 



/// 



For example, uncorrelated tunneling of the two bosons 
would manifest in a transition (1) — > (2) — > (3), i.e., first 
P2 would rise and then P3, while correlated pair tunneling 
would correspond to a direct transition (1) -> (3), i.e., 
P2 would remain zero during the time evolution. Due to 
the normalization of '0(x, y, t), the Pi sum up to unity. 

These truly two-particle quantities are complemented 
by the experimentally easily measurable [37] number of 
bosons in the left well. 



Nt{t) 



dx Pred{x,X,t) = Piit) + ^P2it) , (26) 



which is normalized to a maximal value of one (i.e., Ng = 
1 corresponds to two bosons in the left well) and where 
Predjx, X, t) is the diagonal part of the reduced density 
matrix (20 1. With this normalization, Ng{t) is the two- 



particle counterpart of the single-particle quantity Peit) 
[see Eq. Further measures based on the reduced 

single-particle density matrix are introduced in Sec. |VIB| 



V. DECAY OF TWO INTERACTING 
PARTICLES I: SPECTRAL PROPERTIES 

We now turn to the main object of our present study, 
the tunneling decay of two interacting bosons. We in 
particular wish to clarify which continuum states actually 
support the decay, and under which conditions correlated 
tunneling of a boson pair is observed, rather than inde- 
pendent tunneling of the particles. Taking the symmet- 
ric double- well as starting point, we gradually increase 
the width r of the right well and monitor the transition 



These numerical results perfectly agree with the analytical result 



As is apparent from Fig^a), region (2) consists of two sub- 
regions. Justified by the bosonic symmetry of the wave function 
tpi^, y, t), we integrate solely over one of them, hence in Eq. \2i\ , 
the factor of two appears. 
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FIG. 5. (Color online) Two-particle case: The participation 
ratio PR [magenta, see Eq. ([8|] of the initially left localized 
state IV'(O)) (|17[), as well as its weighted versions (see begin- 
ning of Sec. VA for definition) PRi (black), PR2 (red), and 
PR3 (blue), for diiTerent widths r of the right well, U = 0.2, 
I = 51, 6 = 2, and Vb = 0.1. The dashed cyan curve represents 
a fit with a second-order polynomial y = ax^ + bx + c, see the 
discussion of Eq. | |C3[ ). Note that the width r is sampled at 
larger intervals for r > 1000. The inset shows a magnifica- 
tion for r < 150. At these small values of r, one can clearly 
distinguish peaks of diff'erent width in the PR, which we re- 
fer to as resonances. Since broad resonances also appear in 
PR2, we associate these with single-particle tunneling. Nar- 
row resonances dominate the PR3 and are thus identified with 
two-particle tunneling, see text. As the width r increases, 
the resonances overlap, which renders exclusive one- or two- 
particle tunneling unlikely. 



from tunneling oscillations to the regime of particle loss 
(r ^ £), much as in the single-particle case studied in 
Sec. |III[ In what follows, we will first analyze character- 
istic spectral properties and then, in Sec. |VI[ turn to the 
ensuing dynamics for vanishing and repulsive interaction 
strengths U. 



A. Participation ratio 



In complete analogy with Sec |III[ we analyze the par- 
ticipation ratio PR (and its weighted versions PRi) of 
the initial state. Definitions ([s]) of the PR and ( |Io| ) of 
the PRi directly carry over to the two-particle case. The 
only difference being that, in Eq. ( [lO| , the index is not 
i — r,l but i = 1,2,3, since the configuration space is 
no longer partitioned into left/right, but consists of the 
three distinct regions (1), (2), and (3), see Fig^). Ac- 
cordingly, we replace in Eq. (10 1 the probabihties Pr/ 
with Pi, 2, 3 [see Eqs. ( 23p5 1] ana the {c„} now represent 



the expansion coefficients of the initial state |V'(0)) in the 
two-particle energy eigenbasis {|£'„)} of the double- well 
potential, obtained from diagonalization of Hamiltonian 



([!]). As argued in our discussion of Eq. (10), the PRi 
yield the number of eigenstates that mediate the time 
evolution, weighted with their overlap with the respec- 
tive regions (1), (2), or (3). 

It will prove beneficiary to identify the components of 
the eigenfunctions in terms of the eigenstates of the 
uncoupled double well [Vq — 00). To ease the discussion 
we label the latter as 



\ai, bj 



(27) 



In this notation, a (6) bosons occupy the i-th (j-th) eigen- 
state of the isolated left (right) well, at given interaction 
strength U . The initial state, for example, corresponds to 
|'(/'(t = 0)) — |2i,0), where the component |2i) is given 
by the two-particle ground state ( [T7| ). It is important 
to realize that in a (un-)coupled double-well potential, 
only states of the types |2i,0) and |0, 2^) — where both 
bosons are in one well — depend on the interaction J7, 
while the energy of a state |li,lj) is insensitive to the 
short-ranged contact interaction. Hence, the energy of 
the latter is simply the sum of the single-particle ener- 
gies ef^^ + t-^^^y^A&ad each of the components |li) cor- 
responds to the single-particle wave function \xl) (|xp)) 
see Sec. 1117X1 

All four participation ratios PR and PRi are plotted in 
Fig. [5] versus the width r of the right well, for the initially 
left-localized state IV'(O)) (see Sec. IV) and an exemplary 
interaction strength U = 0.2p^ 

To familiarize ourselves with the quantities at hand, 
we first consider intermediate widths r < 150 of the right 
well (see inset), where one can resolve single resonances in 
the PR. As in the single-particle case, they indicate that 
the initial state |V'(0)) is in resonance with an eigenstate 
of the uncoupled system. Accordingly, we expect pro- 
nounced tunneling only for the resonant configurations 
(see discussion in Sec. IIIB). 

Yet, the appearance of markedly distinct resonance 
widths is striking: Broad resonances occur in PR at 
r « 65,95,130,... and are accompanied by broad max- 
ima of PR2 . In turn, the narrow resonances in PR, found 
at r w 51,75,90,98,113..., go along with marked peaks 
in PR3. The weighted participation ratios unequivocally 
show that broad resonances result from the initial state 
|?/'(0)) being resonant with a state of type |li,li) [re- 
gion(2)] while narrow resonances indicate that the reso- 
nance condition with a state of type |0,2i) [region(3)] is 
met. 

The emergence of different widths constitutes a major 
observation and a fundamental difference with respect to 



12 Here, e^.""' 



is given by Eq. jsjl, but with i replaced by r. 
As long as the interaction is repulsive U > 0, one merely observes 
quantitative differences in the PR, concerning, e.g., the position 
and height of its peaks. For the special case of non-interacting 
particles {U = 0), the PR shows equally spaced resonances of 
similar width as in the single-particle case, but now the three 
states |2i,0), |0, 2^} and participate. 
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the single-particle case (compare Fig. [2]), in which a se- 
quence of equally spaced resonances of monotonically in- 
creasing width was observed (see also Footnote[6]). There, 
we had to consider only resonances between eigenstates 
of the isolated left and right well. In the notation of 
(271 this amounts to resonances between |li, 0) and |0, 1^) 
states. In the two-particle case, not only the resonances 
between |2i,0) and |0, 2^) but also between |2i,0) and 
|li, Ij) states have to be taken into account. 

Given that, the origin of the different widths is ex- 
plained by a simple argument: The resonance width re- 
flects the coupling strength between the participating 
states. Tunneling of a single particle is a first-order pro- 
cess of large coupling strength (and broad resonances), 
while the (two-particle) pair-wise tunneling observed for 
repulsive interactions represents a second-order single- 
particle process of substantially smaller coupling and 
— correspondingly — narrow resonances. [3 More pre- 
cisely, we have found that the two states |2i, 0) and |0, 2^) 
are foremost coupled via off-resonant states of the type 
|li,li). This is corroborated by the observation that, 
for configurations for which the PR displays narrow res- 
onances, also PR2 assumes values in the order of one, 
i.e., also states of the type |li,li) are involved in the 
time evolution. In other words, the observed pair-wise 
tunneling does not correspond to a direct tunneling pro- 
cess from region (1) to region (3), but rather involves an 
intermediate visit of region (2) , although the correspond- 
ing probability P2{t) is small for all times. A similar 
situation is found in the so-called chaos-assisted tunnel- 
ing [58], where tunneling between regular parts of the 
underlying phase space is mediated by quantum states 
which populate the chaotic part of phase space. 

That is, we can tell by the width of the resonance 
whether single-particle or pair-wise tunneling prevails, 
and we confirmed this prediction by the direct time- 
evolution of the initial state, for various configurations. 
The reader is warned though that, in order to observe 
the (second-order) pair-wise tunneling, states of the type 
|li,lj) should be far off-resonant with the initial state. 
Otherwise, the dynamics is dominated by the (first-order) 
single-particle tunneling which, due to the larger cou- 
pling, is a faster process. 

For large widths r > 150, the overlap between consec- 
utive resonances grows and the peaks smear out. Among 
the weighted quantities, PR3 exhibits the strongest in- 
crease and quickly approaches PR, while PRi drops — 
much as in the single-particle case. PR2 increases on 
average and, most importantly, does not drop below a 



^'^ This is consistent with the observation that the pair-wise tun- 
neling (in the symmetric double well) displays much larger oscil- 
lation periods (determined by the coupling strength between the 
wells) than the corresponding single-particle tunneling processes 

For the symmetric double-well, we shall report on this in detail 
elsewhere. 



level of one. This is a first indication that pure two- 
particle loss is unlikely to be observed in the dynamical 
evolution, since the contributing eigenstates have non- 
vanishing overlap with all three regions of configuration 
space. That is, \tp{0)) generally is in resonance with 
|0,2i)- and |li,li)-type states. A priori, this does not 
exclude simultaneous one- and two-particle tunneling. 
Given our above conjecture that single-particle coupling 
is much stronger than two-particle coupling, the former is 
likely to dominate the tunneling process. In the next sub- 
section we elaborate on this point, and further develop 
our spectral picture of the decay of interacting bosons. 



B. Density of states (DOS) 

Let us now attempt to quantify the density of those 
states (DOS) in the quasi-continuum which respectively 
support single- and two-particle tunneling. This is essen- 
tial if we want to convert the participation ratio into a 



decay rate a la Eq. (13), derived above in Sec. IIIB for 



the single-particle case. At the same time, the DOS is of 
fundamental interest since it represents a key ingredient 
in the physics of decaying systems, see, e.g., Ref. [5^. To 
this end, we elaborate on the projections Pi{\En)) of the 
eigenstates [obtained by inserting the |i?„) in Eqs. ( |23[ 
25 1] on the three regions of configuration space as defined 
by Fig. [T|d). For the exemplary case of U — 0.2 (used 
before) and a broad right well r — 3000 {£ = 51), the 
corresponding quantities Pi{\En)) are plotted versus the 
total energy £"„ , in Fig. |6^) . 

Up to a threshold energy oi E = 0.0033, we find that 
PsilEn)) = 1- As the probabihties Pi{\En)) sum up to 
unity, the corresponding eigenstates are entirely localized 
in region (3), i.e., they are of type \0,2i). Above this 
energy, P2{\En)) > indicates eigenstates |li,lj) that 
(partially) populate region (2) where one boson is in each 
well. To zeroth ordeiFj (and in good agreement with the 
numerical results), this threshold energy is given by the 

single-particle ground-state energy in the left well e[^^^ = 
3.79 X 10~^ as given by ([5|. By the same argument, 
the first peaks in Pi{\En)) — corresponding to eigenstates 
partially localized in region (1), i.e., with both bosons in 
the left well — are only observed at about the energy 
of the two-particle ground state in the isolated left well 
ei = 0.0139 [see lowest curve in Fig. [IJa)]. 

We cast this information into a density of states p{E) 
of the environment, i.e., of the right well. Later we will 
rather consider the integrated DOS 



l-E 

n{E) = / dEp{E) 
Jeo 



(28) 



In this estimation of the lowest total energy of a state, 
we neglected the tunneling coupling and the ground-state single- 
particle energy of the right well, which — due to the large width — 



very small compared to e^^''-' , cf. 



Eq. jsjl. 
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FIG. 6. (Color online) Top: Spatial projections Pi{\En)) 
[see Eqs. ( 23 1-(|25|] of the eigenstates \E„) of the asymmet- 
ric double well versus the corresponding energies En, for 
r = 3000 and otherwise same parameters as in Fig. [5] Up 
to approximately the single-particle ground-state energy of 
the left well, ej"*'' — 0.0038 (vertical dashed line), all particles 
are dominantly localized in the right well [region(3)], since 
PsilEn)) = 1. Above this energy, P2{\E„)) > indicates 
eigenstates which (partially) populate region (2), with one bo- 
son in each well. First peaks in Pi{\En)) — corresponding to 
eigenstates with a doubly occupied left well — appear at about 
the energy of the two-particle ground state in the isolated left 
well ei = 0.0139, see Fig.jljs). Middle: The integrated single- 
(red) and two-particle (blue) density of states, for the same 
parameters as in a). Please mind the two y-axis scales. Solid 
lines represent the numerical results ( 29 1 and ( 30 1 , while sym- 



bols correspond to the theoretical predictions (CI I and (C3l 



We note that, for our specific choice of parameters, the peak 
in Pi of panel a) (given by the two-particle energy ei) acciden- 
tally coincides with the kink in n'"^' which results from the 
first excited single-particle state with energy tj^''' = 0.0152. 
Bottom: Only a fraction of the expansion coefficients |c„p of 
the initial state in the double-well's two-particle eigenstates 
(see Eq. Q and Footnote [T]) takes non-zero values in the rel- 
evant energy window around E = ei — 0.0135. The dashed 
line represents a Lorentzian curve of width 7 = 5.9 x 10"** 
iSll. See text for details on the extraction of 7. 



which is less fluctuating. Specifically, we are interested 
in the single-particle quasi-continuum — where only one 
boson escaped from the left well — and in the two-particle 
quasi-continuum — where both particles are located in 
the broad right well. 



The theoretical expressions for the single- p^^iE) and 
two-particle density of states p%^{E) are derived in ap- 
pendix [Cj based on a non-interacting environment. This 
is, of course, a conceptually fundamental assumption to 
be tested in the following. It seems plausible, however, 
to assume that the extremely short-ranged, repulsive in- 
teraction in Hamiltonian ([T]) plays, if at all, a secondary 
role for the energy of the spatially extended continuum 
states. 

In Fi g. [6p ), we compare the theoretical predictions 
(CI), \Ci\ to the numerical results n^^^{E) and 
"-rmm(^)- The latter are obtained from the integration of 
the corresponding curves in panel a) of the same figure, 
explicitly 

<lm{E) - [P2{\EME - E^)] , (29) 
n'LmiE) = m\E.ME - E^)\ . (30) 

m 

For the single-particle quasi-continuum, the agreement 
between theoretical (CI) and numerical results (29 1 is 
very good. The predicted energy gap and consecu- 
tive kinks are clearly observed, merely their position is 
slightly shifted to smaller energies with respect to the 
theoretical values given by Eq. ([s]). This shift results 
from the (tunnel-) coupling b etwe en the wells which is 
neglected in the derivation of ( |Cl[ )p^ 

For the two-particle quasi-continuum, we find almost 
perfect agreement between n^-^^^iE) and n^yij^), al- 
though the latter was derived for non-interacting par- 
ticles. For very large values of the interaction U — \ 
— when the ground state is close to being fermionized — 
we observe that the integrated DOS rt^num is shifted to 
smaller values with respect to n^y^(E^ and thus a gap at 
-E = opens. This shift originates from the interactions 
in the right well [see Sec. IV and Fig.[4]D)] and vanishes for 
r — 00. More importantly, the slope p*^{E) = dn*''/d£', 
remains the same, i.e., the two-particle density of states 
in the quasi-continuum is unaffected by U. Furthermore, 
the gap decreases with increasing width r. This consol- 



idates our initial assumption that a contact interaction 
is of minor importance in the DOS of a very broad well 
which, ideally, stretches out to infinity. 

According to our strategy which proved successful in 
the single-particle case of Sec. |III C[ the DOS together 
with the participation ratios PR and FR^ should yield 
the dynamical decay rate 7 of, say, the particle number 
Ni(t). Inspection of the wave function's expansion coef- 
ficients |c„|^ in the two-particle eigenstates of the dou- 
ble well [see Eq. ([6])] in Fig. [6];) unravels a fundamen- 
tal difference with respect to the single-particle case. In 
contrast to the latter case (see inset of Fig. [2]), only a 



The observed shift diminishes as we decrease the coupling be- 
tween left well and quasi-continuum, e.g., by increasing the bar- 
rier height Vq. 
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FIG. 7. (Color online) Sketch of tunneling decay in the single- and two-particle case. Left: A single particle is prepared in the 
ground state of the left well (energy is marked by the bold black line). For finite barrier height, the ground state couples to all 
eigenstates of the right well within the resonance width (red curve) and the magnitude of the coupling elements is indicated 
by the bold red lines. Middle: Two in tera cting particles are prepared in the two-particle ground stat e (r ed) of the isolated left 
well at energy = (ki i)^ + (^2,1)^ (19 1. The energies (fc^ 1)^ corresponding to the momenta k^ i (18 1 are marked in black. 
For finite barrier height, and assuming mainly uncorrelated tunneling of the particles, the two-particle ground state does not 
couple to all two-particle eigenstates of the right well within the resonance. Instead, the overlap (indicated by the bold red 
lines) is appreciable only for those states which are composed of single-particle eigenstates of the right well with single-particle 
momentum components which match those fc^^i of the initial state. Right: The top panel illustrates a 'matching' eigenstate 
(black) of the right well, i.e., a large overlap, with corresponding (single-particle) energies of the right well marked in blue. The 
bottom panel illustrates a 'non-matching' eigenstate of the right well, where the corresponding single-particle energies (blue) 
are highly off-resonant with respect to the energies (fc^^i)^ of the left well (black). 



fraction of the c„ take a finite value, within the relevant 
energy window around the two-particles ground-state en- 
ergy E = d = 0.0139. As a result of the vanishing |c„p 
the two-particle DOS (C3) overestimates the number of 



effectively contributing states. This impedes to proceed 
as in the single-particle case along Eqs. (Ill, (12): To 



identify 7 with the quotient of participation ratio and 
(in this case overestimated) DOS would underestimate 
the decay rate 7. Nonetheless, we extract a resonance 
width 



— unlike the single-particle case — it is not enough to con- 
sider the total energy e„ of the initial (two-particle) state 
alone, as the greater part of energetically matching two- 
particle states of the right well actually do not couple to 
the initial state, due to the mismatch between the under- 
lying single momenta. Put differently, tunneling decay in 
the two-particle case is not about resonances between the 
total two-particle energies but, rather requires the indi- 
vidual matching of the two momenta fc^^i in the left well 
to the momenta in the quasi-continuum. 



5.9 X 10" 



(31) 



from the 50% quantiltf^ of the expansion coefficient's 
distribution |c„p, see Fig. [gJ:), and comment on its role 
for the dynamical evolution in the next section. 

What, at first sight, appears to be a further compli- 
cation, de facto much elucidates the intricate tunneling 
mechanism: We have found that the vast majority of van- 
ishing components within the resonance width originate 
from quasi-continuum states that do have the appropri- 
ate two-particle energy of the right well e„ « e\^^^ + ej"^'' 

[^but their constituting single-particle energies e^*^^ do 
not match the momenta kij that enter the energy ei 
([19]) of the initial state |'i/'(0)). This demonstrates that 



Specifically, we evaluated 7 from the condition 
kriprfSn = 0.5. Here, the center of the distribu- 
tion, xq, is determined by /^J^ \cn\'^dE„ = 0.5. 
^® If we ignore the influence of the repulsive interaction in the envi- 
ronment as above, the true two-particle energy e„ can be replaced 

(sp) 



by a sum of two single-particle energies of the right well e 
given by Eq. ||5|, but with £ replaced by r. 



In the latter sense, we are dealing with a case of single- 
particle tunneling, for which the inter-particle correla- 
tions play a minor role. At the same time, this matching 
principle — which we schematically illustrate in Fig. [7] — 
is of genuine two-body nature as it crucially relies on 
the (increased energy) of the two-particle wave vectors. 
To properly interpret the illustration, we note that the 
assumption of uncorrelated single-particle tunneling im- 
plies that the initial |2i, 0)-state is not directly converted 
into two-particle quasi-continuum states, but via inter- 
mediate states of type |li, 1^). 



Let us summarize the main message of the spectral 
analysis: For moderate values of the width r, isolated 
resonances exist in the PR and pair-wise two-particle 
tunneling from the initially prepared ground state |^/'(0)) 
can be expected for appropriate choices of the system 
parameters r, £, and U . In the quasi-continuum case of 
r -> 00 and repulsive interactions, resonances strongly 
overlap, compare Fig. [5j Although only a fraction of 
the two-particle continuum states was found to actually 
support the particle decay, the growing PR2 and PR3 in- 
dicate that — irrespective of r — the initial state couples 
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to both, the single- and the two-particle continuum p"] 

We therefore develop the intuition that states which 
are transformed into each other by one single-particle 
process are considerably stronger coupled than those 
which require a second-order process. In total, we thus 
expect (second-order) pair-wise loss to be less generic 
than uncorrelated first-order tunneling to the quasi- 
continuuni. 



VI. PARTICLE LOSS OF TWO INTERACTING 
BOSONS II: DYNAMICS 

Up to this point we have not discussed the time evolu- 
tion of two-particle loss, a topic which has only recently 
gained momentum. First experiments range from atom 
loss off large BECs, induced by electron scattering [5l] 
or photo ionization |15j to the tunneling-induced loss in 
highly controlled few-fermion systems [37 . Only a hand- 
ful theoretical works treat the full-fledged microscopic de- 
cay problem considered here, and they rely on the propa- 
gation of the time-dependent Schrodinger equation either 
by direct integration ,39 , by combining a matrix-product 
state approach with Bose-Hubbard-like chains [ID] , or via 
multi-configurational Hartree-Fock methods [U . 

In this section, we complement the results of our spec- 
tral analysis by the associated time evolution which we 
obtain by the spectral decomposition ([7| of the initial 
state ( 17) [but now for the two-particle quantities, result- 
ing from the exact diagonalization of Q], and comment 
on existing results as we proceed. We use the same ini- 
tial state as in the spectral analysis (see Sec. IV B I, and 
the same width r = 3000 of the right well, which we also 
considered in our analysis of the single-particle dynamics 
in Sec. Ull 



A. Decay of the particle number Ne{t) and the 
two-particle probabilities Pi{t) 

We start the discussion with an experimentally easily 
measurable quantity: the normalized number of bosons 



Ni{t) in the left potential well (26). This quantity is plot- 
ted in Fig. [8^), on a semi- logarithmic scale, for various 
values of the inter-particle interaction strength U. For 
U = 0, Ng{t) decays exponentially with the decay rate 



of the single-particle ground state 7{*^'' ( |15[ ) (indicated 
by the dotted line). This is not surprising since, in the 
interaction-free case, the system should effectively reduce 
to a single-particle problem]^ As U increases, however. 



20 -^Yg remark that this is in perfect agreement with Fig. |6p) which 
shows that, at the energy ei = 0.0139 of the initial two-particle 
ground state in the left well, |0, 2^)- and |li, li)-type states are 
populated in the continuum, i.c, both are in resonance with 

\m)- 

Effects that arise from the indistinguishability of (non- and in- 
finitely strong interacting) particles are discussed, e.g., in |60H63| . 



pronounced deviations from a (straight-line) exponential 
behavior arise, especially for times t < 5000. The second 
striking feature is that all curves asymptotically assume 
the same slope. 

Insight is provided by the two-particle probabilities 
Pi=i,2,3{t) = Pi{\ip{t))) which — by virtue of the defi- 
nition (26) — determine Ng{t). For the sake of visual 



clarity, we present the results for the exemplary value of 
U ^ 0.2 in Fig. Initially, the probability Pi(t), to 

find both bosons in the left well, drops. At the same 
time, the probability P2{t) that one particle is located 
in the left and one in the right well increases. Only 
after that does the probability to find both bosons in 
the right well, Psit), rise. The key observation 
is that P2{t) becomes considerably large, independently 
of the interaction strength U. We conclude that, to 
a large extent, the bosons do not leave the trap as a 
pair — corresponding to a direct transition between re- 
gions (1) — )■ (3) — but sequentially leave the trap, i.e., 
through regions (1) — > (2) — > (3). This confirms our pre- 
diction based on the spectral analysis of the preceding 
section which we further substantiate in the following: 
Namely, for U > 0, the predominant loss mechanism is 
uncorrelated, i.e. independent, tunneling of the bosons 
[321 mi m] albeit the presence of interactions. 

With the help of the Pi (t) , we can now readily explain 
the short- and long-time behavior of the particle num- 
ber Ni{t). The former will be essentially determined by 
the probability Pi{t) which (in agreement with [33]) is 
found to decay exponentially for all values of U. We re- 
mind the reader that the initial state of the two bosons 
'0fci,ifc2 1 (2;, y) is characterized by the two wave vectors 
fci^i and ^2,1- Encouraged by the result of our spectral 
analysis and the behavior of the Pi{t), we assume an in- 
dependent particle picture as in [3^] , and associate a sep- 
arate decay constant 7^. ^ with each wave vector. That 
is, the probability Pi{t) to find both particles in the left 
well is given by the product 



Piit) 



(32) 



where the ki^i depend on U via Eq. (18) 



Let us first consider the trivial case U = Q. Here, 
both particles are in the single-particle ground state 
and the decay is determined by the corresponding single- 



7fc2 



7i 



(sp) 



particle rate given by (151, i.e., 7fc^ 
For intermediate values 0< f7 < 00, the iJ-dependence 
of the rates can be inferred from Fig. l4k): The momen- 



tum ^2,1 is largely independent of U and given by the 
single-particle ground state momentum, hence we have 
7fc2 1 ~ 7i*^^- In contrast, the momentum fci i — and 
thus the associated energy k1 ^ — grows with U , leading 
to an increase of 7^^ ^ (since the effective barrier height 
seen by the particle is decreased). This growth continues 
until, at large values of J7, the system is fermionized and 
7fc^ J takes the value of the decay rate 72*^' associated 



with the single-particle first excited state (16). In total. 



the decay of Pi{t) therefore increases with U and takes 
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FIG. 8. Particle loss of two interacting bosons in an asymmetric double-well potential, in which the broad right well emulates 
the continuum. The time evolution is obtained from diagonalization of ([T]) with / — 51, b — 2, and r — 3000, and the particles 
are initially prepared in the ground state of the isolated left well (see Sec. |IV[ ) Left: The normalized particle number Ne{t) 
(26 1 in the left well, versus time t, for various values of the inter-particle interaction strength U. For vanishing interaction, 
an exponential decay is observed. As the interaction increases, pronounced deviations from the interaction-free exponential 
behavior arise. To guide the eye, the dash-dotted (dashed) line indicates an exp one nt ial d ecay with a rate corresponding to the 
single-particle ground (single-particle ground plus first excited) state [see Eqs. (15 I, (161, and (33l]. Right: For the exemplary 
interaction strength U — 0.2, we plot the probabilities Pi{t) — Pi{\ilj{t))) to find the particle in the three regions (l)-(3) of 
configuration space corresponding to two, one, and zero bosons populating the left well, respectively [see Eqs. ( |23| l-(|25[l and 
Fig. [ija)]: The strong initial increase of P2{t) compared to Psit) indicates that the particles predominantly tunnel in the 
sequence (1)— >-(2)— >(3). That is, they tunnel sequentially and not as a correlated boson pair. 



the limiting cases: 



Piit) 



[/ = 0, 
[/ = oo 



(33) 



One thus finds that, for short times t, the (for non- 
vanishing U larger) rate i dominates the time evo- 
lution of Ng{t), and defines an upper limit (7^**^^ +72*^'') 
for the initial decay rate [indicated by the dashed line in 
Fig#)]. ^ 

Concerning the long-time evolution, we assume that, 
after the first boson left the trap, the remaining one 
populates the single-particle ground statej^ Hence, this 
second boson decays with the single-particle decay rate 

7i 

erned by 7^*^^ i.e. 



v(*p) Put differently, the loss out of region (2) is gov- 



P2(t> 1) oc e 



(34) 



In agreement with our numerical simulations, the U- 
independent, smaller decay rate ~ 7^2 ^ (dash- 

dotted line) governs the asymptotic behavior of Ni{t ^ 
1) for all values of f7 > 0. 

From the preceding discussion, we cannot but conclude 
that the assumption of uncorrelated tunneling very well 



This assumption is substantiated further down, in Sec |VI b] in 
the analysis of the reduced density matrix in energy space. 



reproduces the observed decay of the particle number 
inside the left well. 

We point out that our results on Ng{t) and Pi{t) 
excellently agree with those of Ref. in which a 

delta- like potential barrier was considered. In contrast, 
Refs. im [B3] do not mention such non-exponential 
decay of the particle number Ni{t). Neglecting the dif- 
ferent trap geometry, the parameters used in the main 
part of |42j roughly correspond to interatomic interac- 
tions of the order oi U ~ 0.05; inspection of Fig. |4^) 
reveals that, for U sa 0.05, the difference in the wave 
vectors ki i is comparatively small. Thus also the rates 
7fcj^ J and 1 barely differ, making it difficult to ob- 
serve non-exponential decay. In the associated Support- 
ing Information of [42j , the time evolution for two bosons 
and a seven-fold larger interaction strength is reported 
in Fig. S2. A careful re-analysis of the data (semi- 
logarithmic plot) reveals that the decay indeed consists 
of two exponentials, as discussed above, leaving no in- 
consistencies. 

We note in passing that we determined the decay rate 
of Pi{t) by an exponential fit to be 



^U=0.2 ^ g 8 X ^0-4 



(35) 



which is in good agreement with the spectrally extracted 
width (31) of the Lorentzian curve in Fig. l6b). As we 



explain in Appendix [D| this indicates that for our setup, 
the probability Pi(t) corresponds to the survival proba- 
bilityP,urv = |(^(0)|V(t))|'. 
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FIG. 9. (Color online) Reduced single-particle quantities of the two boson decay. Left: Integrated diagonal d{x,t) (361 of the 
reduced single-particle density matrix p^cdi^, x,t) versus position x, for different evolution times t and inter-particle interaction 
strengths U, for otherwise same parameters as in Fig.|8] For vanishing interactions, the particles leave the left well independently 
and with the same probability. Hence, one observes a continuous increase of d{x,t). In contrast, for U — 0.2, the bosons tunnel 
to the right well at different rates and different asymptotic momenta. This leads to a spatial separation in the quasi-continuum, 
manifest in the plateau around the interval x £ [1500,2400] (for t = 15000). Right: The diagonal pred{E, E,t) (37 1 of the 
reduced single-particle density matrix in energy representation versus the energy E, for the same parameters and color code as 
in panels a) and b). In the interaction-free case U — 0, one observes a dominant Lorentzian-shaped peak at (dashed vertical 
line) which comprises about 98% of the total probability and a smaller peak at e'^^^ (dash-dotted vertical line). The reduced 
probability density pred{E, E,t) does not change in time, hence all curves fall on top of each other. For repulsive interactions, 
the initial distribution is qualitatively similar to the case (7 = 0. During intermediate times t < 6000, the interaction energy 
is transformed into kinetic energy and the distribution picd{E, E,t) changes. As a result of this conversion process, two peaks 
emerge, located at the energy associated with the two wave vectors kf i ~ and fc2,i ~ — ej^"''' (dash-double-dotted vertical 
line) which constitute the initial state. Each of the peaks in the asymptotic distribution comprises 50% of the probability, what 
corroborates the picture of two asymptotically independent particles. 



B. Reduced density matrix in coordinate and 
energy space 

We now wrap up the discussion of the decay dynamics 
with two complementing reduced single-particle observ- 
ables. These are potentially easier accessible in the ex- 
periment than e.g., the two-particle probabilities Pi{t). 
On the other hand, they add to the understanding of the 
time-resolved tunneling process. We start with the in- 
tegrated diagonal d(x, t) of the reduced density matrix 
/9red(a;,a;), i.e., 



d{x,t)= i c\x' p^cd{x' ,x' ,t) . 
Jo 



(36) 



which denotes the probability to find — at a given time 
t — a boson in the interval [0, x] of configuration space. 
Consequently, as long as d{x, t) increases with x (for fixed 
t), there is a non- vanishing probability that a boson re- 
sides at the position x. As will become clear in the next 
paragraph, when seen as a function of time, d{x, t) visu- 
alizes how the bosons leave the left well. We finally note 
that d{£,t) — N({t), as the left well stretches from x = 
to X — i. 

In the left panels of Fig. |9j we plot d{x, t) versus the 
position X, for various times t and two representative 
values of U . To familiarize ourselves with the quantity 
at hand, we first consider the interaction-free case J7 = 0, 



shown in panel b). For example, (i(51, 15000) ~ 0.2, that 
is, at time t = 15000, about 80% of the particles have 
tunneled through the potential barrier located at x = 
51. More importantly, d{x > 51, t) shows a monotonous 
increase in x for all (fixed) times. Since for U — 0, the 
two bosons are in the same single-particle (ground) state 
[see Eq. (21)] and do not interact with each other, this 
constitutes our reference case of two bosons leaving the 
left well with the same probability. 

Contrast this with repulsively interacting bosons [U = 
0.2, panel b)]: As time evolves, a plateau at d{x, t) « 1/2 
emerges, see, for example, the magenta curve {t = 15000) 
around x € [1500,2400]. According to the definition of 
d{x, t), this implies that the probability to find a particle 
within this plateau is about zero. Put differently, with 
probability one half one finds a boson between a; = 
and X — 1500 and equally likely a boson can be detected 
between x = 2400 and x = 3000. This behavior nicely 
fits our above developed picture, in which the bosons 
independently leave the trap but (for U > 0) with differ- 
ent loss rates 7^-1. For even larger times, this plateau 
is stretched and the two particles are further separated, 
indicating that they travel at different speeds. 

To corroborate this statement and to complement the 
discussion on the particles' spatial distribution, we finally 
discuss the fingerprint of the particle loss in energy do- 
main — let aside the advantage that the energy (i.e., the 
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momentum) of a single particle may be easier to mea- 
sure than its position. For this purpose, we consider the 
reduced single-particle density matrix in energy represen- 
tation 

p„d(S, E', t) = T{E, x)p,ed(x, y, t)T\E', y) , (37) 

T being the unitary basis transformation from configu- 
ration space to the energy domain. 

In what follows, the reader is warned that the reduced 
quantities are associated with measurements on a single 
particle. While even a single boson carries some infor- 
mation on the total system (due to symmetrization) , the 
reduced quantities do not reflect correlations between the 
particles, e.g., the probability that both particles have 
the same energy, cannot be deduced from picd- Here, 
we focus on the diagonal elements prcd{E,E,t), i.e., the 
probability density to find one particle with energy E. In 
the right panels of Fig. |9]we present pred{E,E,t) versus 
E, for the same parameters as in panels a) and b) of the 
same figure. 

Paralleling the discussion for d{x,t), we take the 
interaction-free case [panel c)] as our starting point. Ini- 
tially, the two bosons populate the left well and the two- 
particle initial state |'0(O)) is characterized by the two 
momenta ki^i and fc2,i. For U — 0, the latter are identi- 
cal and furthermore coincide with the momentum fcj"^^ of 
the single-particle ground state, see Fig. |4]a,). Intuitively, 
one would thus expect to observe a peak in Prod(-E', E, 0) 
at the corresponding energy e^"^' (dashed vertical line). 
Indeed, we identify one dominating, Lorentzian-shaped 
peak (comprising about 98% of the probabilitjj^ , lo- 
cated slightly below e"i^\ As we already remarked in 
Sec. |III[ this energy shift results from the coupling of the 
left well to the quasi-continuum and increases with the 
energy of the initial state [see also Footnote 17 . The 
remaining probability is contained in the second, inferior 
peak, situated close to ej*^'' (dash-dotted vertical line). 
We remark that this peak would as well appear in a truly 
single-particle calculation and results from a finite over- 
lap of the initial state \ip{0)) and quasi-continuum states 
that also involve excited single-particle states of the left 
well. These minor corrections aside, intuition does not 
fail us. More importantly than the position of the peaks, 
Picd{E, E, t) does not change its shape as time evolves: all 
curves exactly fall on top of another. Accordingly, both 
bosons follow the same energy distribution at all times, 
that is, they move with the single-particle ground state 
momentum through the quasi-continuum, as expected for 
the non-interacting case. 

How do repulsive interactions alter this picture \U = 
0.2, panel d)]? In contrast to the case C/ = 0, the ini- 
tial momenta fci_i and ^2,1 no longer coincide with the 



single-particle momenta k!f^^ (i.e., with the e^'*'^'') and 
one is tempted to expect a direct manifestation thereof 
in the initial distribution p^cd{E,E,Q). Yet, the differ- 
ence in pred (-£',£', 0) compared to the case C/ = seems 
marginal: For t = 0, we still observe two peaks at the 
same positions as before, save that the main peak cen- 
tered around e^"^' now comprises only about 90% of the 
probability!^ The reason why the kij do not emerge 
is that, at i = 0, only single-particle eigenstates of the 
left well can be populated on the reduced single-particle 
level. Only after the particles tunnel to the right well, 
quasi-continuously distributed states are available, and 
every energy component can contribute to the dynamics. 
In accordance with the warning issued after Eq. ( 37 1 , 



The probability comprised by each peak is calculated via inte- 
gration of pj-f,^{E,E). As integration boundaries, we take the 
middle between the centers of neighboring peaks. 



Prod(£', E, 0) rather reflects the interaction- induced finite 
contribution of IX2) to the initial state |V'(0)) than the 
(truly two-particle) momenta fci_i. 

The influence of the interaction becomes fully appar- 
ent as time evolves: The peak at e'*^'' declinesr^ while 
the one at £2^^^ entirely vanishes, and a new one centered 
around the energy E — 9.7 x 10"'^ emerges. The change 
in the shape of p^cd{E, E, t), clearly indicates that energy 
is being redistributed among the particles. In the inde- 
pendent decay picture, the particle that remains in the 
left well populates the single-particle ground state, i.e., 
it has energy e^'*''^ w (^2,1)^- Due to energy conservation, 
the boson that tunneled out of the left well should carry 
the energy ei — e[''^^ ~ (/ci^)^. Indeed, we can associate 
the emerging peak with the latter energy, while the first 
peak still corresponds to the ground state energy e[''^\ 
in good agreement with the result of |42| . Our interpre- 
tation that the two peaks represent the two bosons is 
corroborated by the fact that (for times t > 6000) each 
peak comprises half of the probability. We furthermore 
found that the energy widths of the two (asymptotically) 
Lorentzian peaks in pj-^diE^ E, t) approximately equal the 

two decay constants 7^ and 7^"^' associated with the 
decay of the two particles. All three observations con- 
cerning the position, integrated area, and width of the 
peaks in p^f.d{,E,E,t) further substantiate the picture of 
uncorrelated, independent decay of the bosons. Hence, 
PicdiE, E, t) represents an experimentally easily measur- 
able quantity which bears substantial information on the 
tunneling decay. 

Let us now take a closer look at the energy conversion 
which takes place during the tunneling process. Initially, 
both bosons populate the left well, and interact with each 
other. Once they tunnel to the right well, their interac- 
tion energy is converted into kinetic one. This process 



The largest part of the remaining probability is contained in 
the second peak [which is increased with respect to panel c)] 
while higher than first excited single-particle eigenstates con- 
tribute with considerably less weight. 

This decline is less pronounced on the logarithmic scale, but can 
be inferred from the increasing second peak. 
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lasts until t « 6000, until when the second peak grows 
and takes a Gaussian shape [see red curve in Fig. [9]:)] . For 
larger times {t > 6000), one may imagine that the boson 
with the larger decay rate has completely tunneled to the 
right well and is moving at a larger momentum ki i. 
This is also approximately the time when the plateau in 
d{x,t) appears, see panel b). Then the bosons cease to 
interact, and the peak takes its asymptotic Lorentzian 
shape. As a result, the interaction part of the total en- 
ergy ([T]) effectively vanishes at these large times and we 
find 

{i^it)\Hsp 1 + 1(8) Hspli^it)) = 2 ti{HspP,ed) = ei . 

(38) 

That is, the total energy ei, determined by the initial 
state, is given by twice the average single-particle energy 
what is certainly not the case for smaller times, when the 
particles still interact. 

As a last remark, we have numerically confirmed that, 
for r — >■ cx), the interactions between the bosons in the 
right well have almost no effect on the loss dynamics 
Hence, it is well justified to ignore them, what is in per- 
fect agreement with our results from Sec. |VB[ where 
both, the single-particle and the two-particle density of 
states of the quasi-continuum proved almost independent 
of U. While we stress that this result was derived only 
for 5-like interactions, the possibility to neglect interac- 
tions in the environment greatly simplifies perturbative 
approaches to the loss dynamics, such as a master equa- 
tion ansatz. 

Summarizing this subsection, we have found that re- 
pulsive interactions strongly modify the loss dynamics 
and lead to a multi-exponential decay of the boson num- 
ber. Notwithstanding, the dominant decay mechanism 
is found to be uncorrelated tunneling of single parti- 
cles. That is, in none of the studied examples, signifi- 
cant correlated two-particle loss was observed. Recalling 
the spectral analysis of the previous section, both single 
and two-particle tunneling were found to be energetically 
possible. Yet, the simple physical picture developed in 
Sec. |V] turns out to be correct: single-particle tunneling 
— as a first order process — wins over the slower two- 
particle process of second order and thus constitutes the 
dominant tunneling mechanism. 

VII. SUMMARY AND DISCUSSION 

In the present work we studied the tunneling of two 
ultracold bosons, initially prepared in the left site of 
a double-well potential, as the width of the right well 
was gradually increased in subsequent realizations of the 
trap geometry. Seen from the left well, a broad right 



We set f/ = in the right well, and found almost no difference 
with respect to the curves of Figlsl 



well mimics an unconfined configuration space (i.e., an 
environment) with a dense quasi-continuum of states to 
which the particles escape. Through the numerically ex- 
act diagonalization of the full two-body Hamiltonian, we 
could identify those quasi-continuum states which actu- 
ally support this decay and expose the role played by 
the repulsive interactions. We could also conclude that 
the many-body decay process is governed by indepen- 
dent tunneling of the bosons rather than by tunneling of 
a boson pair. We briefly review the main points. 

As a starting point and frame of reference, we first 
investigated the single boson case. In the participation 
ratio PR of the initial state, we found a simple spectral 
tool, to predict particle dynamics: Its maxima with re- 
spect to the width r of the right well, termed resonances, 
reliably indicate those trap geometries for which substan- 
tial tunneling occurs. Still isolated at comparable widths 
of the two wells, the resonances increasingly overlap as 
the quasi-continuum r — >■ c» is approached. This was 
identified as a spectral signature of the transition from 
tunneling oscillations to tunneling decay. The rate of 
the corresponding exponential particle loss was extracted 
from the asymptotic growth of the PR and found to be 
in good agreement with the associated time evolution as 
well as with analytical results. 

In the main part of this work, we then analyzed the 
case of two bosons, initially prepared in the two-body 
ground state of the isolated left well. As a first ma- 
jor difference compared to the single-particle case, the 
participation ratio PR exhibits resonances of two differ- 
ent widths. We introduced the concept of configuration- 
space sensitive participation ratios PRj and unequivo- 
cally identified narrow resonances with second-order pair- 
wise tunneling, and broad resonances with uncorrelated 
first-order tunneling. 

In the quasi-continuum limit of a broad right well, the 
resonances of each kind were found to overlap, what im- 
plies that a priori both, single-particle and two-particle 
states of the continuum are available for the loss dynam- 
ics. This statement was corroborated by the excellent 
agreement of our analytically and numerically evaluated 
single- and two-particle density of states which further- 
more proved to be independent of the inter-particle in- 
teraction U. 

As a second major difference with respect to the single- 
particle decay, we showed that not all energetically al- 
lowed quasi-continuum states contribute to the tunnel- 
ing, but only a small fraction thereof. It was found 
that this fraction is determined by the condition that 
the single-particle energies, which underlie the quasi- 
continuum, match the energies associated with the two- 
particle momenta fc^^i of the initial state. 

We complemented our spectral analysis with the study 
of the associated time evolution. In agreement with pre- 
vious works [39], we found that depending on U, the 
initial decay of the boson number in the left well Ng{t) 
fundamentally deviates from the (interaction-free) ex- 
ponential behavior, while the asymptotic decay of two 
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bosons was found to be exponential and interaction in- 
dependent. The results were explained in terms of the 
interaction-dependent two-particle wave vectors /c^ i. We 
furthermore showed that the essential characteristics of 
the two-particle loss process can be extracted from mea- 
surements of reduced single-particle quantities and elab- 
orated on the conversion from potential to kinetic energy 
during the tunneling process. 

In combining the results of the static (spectral) and dy- 
namic analysis, we developed a simple and clear picture 
of the tunneling decay: The repulsive interactions first 
and foremost modify the energy of the bosons what even- 
tually determines their tunneling rate and final momenta. 
That is, although the initial state is of truly two-body 
nature, the tunneling decay is well reproduced by associ- 
ating the two-particle wave vectors fc^^i with independent 
particles. Furthermore, through appropriate choice of the 
system parameters, (second-order) tunneling of a boson 
pair may be enforced for moderate asymmetries of the 
trap where the resonances are still isolated. The particle 
loss, however, is dominated by single-particle tunneling 
which (as a first-order process) wins over the slower two- 
particle process of pair-tunneling. 

Experimentally, quantum systems of few ultracold par- 
ticles have recently been realized for arbitrary values of 
the inter-particle interactions [TT] and the loss process of 
the first particle has already been monitored ^37^ . Al- 
though even box shaped potentials seem experimentally 
feasible [8l-fT0|. our observations should not qualitatively 
depend on the specific form of the trapping geometry. 
Hence, an experimental investigation of the few-body de- 
cay is realizable with state-of-the-art technologies. Be- 
sides the decay, the isolated resonances in the participa- 
tion ratio for comparable widths of the two wells demon- 
strate an experimentally feasible way to address — in a 
controlled way — two- and single-particle tunneling by 
changing the geometry of the trap, rather than tuning 
additional magnetic fields. 

In closing, let us outline future directions of this 
research area. Firstly, the fraction of pair-tunneling 
in the loss dynamics — as low as it may be — should 
be determined in a future study together with a pre- 
scription to experimentally detect it. Secondly, due to 
their inherent symmetry, bosonic systems exhibit particle 
bunching j6Sj . even for massive, non-interacting bosons 



[551 [57] . One may thus very well ask to which extent 
the interaction-induced fermionization process observed 
in the present setup modifies the bunching behavior. In- 
timately related is the dependence of the dynamics on 
initially prepared states different than the ground state 
of the left well, e.g. excited states or non-eigenstates 
pOl [55] . Thirdly, we here considered a 5-like interac- 
tion which represents the simplest two-body interaction. 
What will happen when more complex (e.g., long-range. 
Coulomb) interactions come into play? First results exist 
for the so-called escape dynamics j69j (a precursor of tun- 
neling decay in which the confining potential is suddenly 
switched off and particles escape to free space) . A recent 
study [70] indicates that charged ultracold bosons in con- 
tinuous potentials may exhibit quantum chaotic behavior 
(much as the one observed in quantum lattice models, 
see, e.g. [32]), which may alter as well the decay process 
[71] . Similarly, the study of attractively bound boson- 
pairs [75] would be of both experimental and theoretical 
interest. E.g., will the tunneling break the bonding or 
will the particles tunnel as a pair? 

We finally point out a long-term objective of eminent 
computational relevance. While the decay of a small 
number of bosons can still be treated exactly [12j or 
within a Bose-Hubbard-type modeling |40] . the micro- 
scopic dynamics of larger boson numbers is beyond nu- 
merical reach. The observed (predominantly) uncorre- 
lated single-particle tunneling decay might advance the 
development of perturbative approaches, such as a mas- 
ter equation. Without resorting to heuristically intro- 
duced approximations, one should base the master equa- 
tion ansatz on the energy-dependent loss rates described 
above and compare the outcome to existing results. 
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Appendix A: Two-body Hamiltonian 

The most general two-boson Hamiltonian reads: 
H ^ j dx ■^\x)Hsp-^{x) 

dx I dy'^\x)'^'^{y)W{x,y)-^{x)'^{y) , (Al) 



where '^'^^\x) are bosonic field operators which anni- 
hilate (create) a boson at position x and Hsp is the 
single-particle Hamiltonian defined in Eq.([2|. The sec- 
ond term accounts for the interaction W{x^y) between 
the two particles with coordinates x and y, respectively. 
We construct the two-particle Hamiltonian by expand- 



ing the bosonic field operators in (Al) in terms of the 
eigenfunctions E')j^^\x) = (xji^j"^^) of the single-particle 



Hamiltonian Hsp, i.e 
*(x) 



Y,a,E^;^\x). 





(A2) 



This yields 



where E'j''^'' are the single-particle eigenenergies and gij 
is the matrix element that accounts for the contact inter- 
action W{x,y) = U5{x — y) (see discussion after ([3| in 
the main text) between the particles, i.e., 

9^jim = ul dx [El'^\x)E^;^\x)y El^^\x)EilPH^) ■ 

(A4) 

In our numerical treatment, we retain the lowest ricut 
single-particle eigenfunctions leading to a two-particle 
Hilbert space of dimension {ncut){ncut + l)/2 which 
is spanned by the basis {S{\Ei7'^) (g) jE'i"^^))}. Here 
{\Em^^)} is the single-particle energy eigenbasis, S is the 
symmetrization operator, and n < m = l,...,ncut- In 
this approach, convergence is reached when the results 
do not change upon further increase of Ucut- Typically, 
several hundred single-particle states are taken into ac- 
count which enables us to study a broad range of repul- 
sive interactions U. We note that this treatment goes 
far beyond the single-band approximation often assumed 
in the context of ultracold atoms, e.g. within the Bose- 
Hubbard approach [73]. 



"•j^j + E 9tjima\a\aiam , (A3) 



Im 



Appendix B: Analytical treatment of single-particle 

loss 



For an infinitely broad right well (r — cxj), we can 
readily formulate the single-particle loss as a scattering 
problem. In each of the three intervals I, II, and III of 
configuration space [see Fig. [T^)], we make a (textbook- 
like) plane-wave ansatz for the single-particle wave func- 
tion: 

{Aexp(— ifcx) -|- B exp{ikx) x Gl 
C exp{— Kx) + D ex.p{Kx) xell , (Bl) 
F exp{—ikx) + G exp{ikx) x £ HI 

with = Vo — k'^. The postulated continuity of ip{x) and 
its first derivative at the borders of the intervals yields 
four equations for the six coefficients A to F. We are 
seeking purely outgoing solutions and thus fix the ampli- 
tude F{k) = 0. The latter equation has solutions kn^^ 
that lie in the complex plane, leading to complex energies 
En^'' = {kn^'^)'^ whose imaginary part corresponds to the 
decay rate of the n-th eigenstate \Xn) of the isolated left 
well. 



Appendix C: Two- and single-particle density of 
states 

In this appendix, we derive theoretical expressions for 
the single- and two-particle density of states in the quasi- 
continuum, what corresponds to one and two particles 
having escaped from the left well, respectively. The com- 
parison to our numerical data is found in Sec. |VB[ see 
discussion of Fig. [6j 

We assume that the single-particle quasi-continuum is 
formed by product states \Xn)^\Xm) of the single-particle 
eigenstates of the left and right well. This assumption 
neglects the tunneling coupling and is thus strictly valid 
only for uncoupled wells. The interaction between bosons 
in the left and right well, on the other hand, can be safely 
ignored since we are dealing with a contact potential in 
Hamiltonian (jlj. According to Eq. (12), the |Xm) deter- 
mine the integrated single-particle DOS of the right well. 
Since we would like to obtain the integrated DOS as a 
function of the two-particle energy E, we have to addi- 
tionally take into account the energy of the states jxl)- 
By summing over all combinations {n, m} with total en- 
ergy below a fixed value E, we obtain the integrated DOS 
of the single-particle quasi-continuum, 



E 



E - eif ) e{E - 



Asp)' 



(CI) 



where the em^^ are the eigenenergies of the isolated left 
well [see prior to Eq. ^] and 6{-) is the Heaviside step 
function. We note that the latter implies an energy gap 
of size e^'^'' and potentially also discontinuities (kinks) at 
consecutive energies en^"^ ■ 
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For the two-particle continuum, an analytical expres- 
sion is more involved, as the two-particle eigenenergies 
result from ^l^ and the coupled equations (flSl). In- 



stead, let us assume for the nonce that the contact in- 
teraction plays a secondary role as far as the energy is 
concerned, since the right well consists of the spatially 
extended quasi-continuum states. We shall discuss the 
validity of this assumption in the main text. Then, for 
two non-interacting bosons, the associated two-particle 
density of states is given by [HI [75] : 



(C2) 



where ★ indicates the convolution and p^'*^^ is the single- 



particle density of the right well, see definition ( 12 ) with 
i replaced by r. After a short calculation, we obtain 
the integrated two-particle density of states in the quasi- 
continuum, 



An 



-E 



27r 



/2E 



(C3) 



Two remarks on this result are in order: Firstly, 
Eq. (C2| predicts a second-order polynomial increase of 
p%{E) with r. Recalling the behavior of PR in the quasi- 
continuum limit of the single-particle case [see around 
Eq.([l3|], we deduce that the PR in the two-particle prob- 
lem should grow with the same functional dependence on 
r. This is confirmed by the fit with a second-order poly- 
nomial in Fig. [5] Thus, values of r > 1000 should be well 
in the quasi-continuum, i.e., suitable to study particle 
loss. Secondly, we note in passing that in the limit of an 
infinitely large well r — > 00, one finds n*^(_B) oc E. Not 
surprisingly, the resulting constant DOS for two (non- 
interacting) particles in a one-dimensional system is ex- 
actly the DOS of one particle in two dimensions. 



Appendix D: Weak coupling regime and survival 
probability 

We here discuss the link between the probability to 
find all particles in the left well and the so-called survival 
probability 



,{t)^\{m\mw 



(Dl) 



which denotes the probability to find the system at a 
given time t in the initial state The former prob- 

ability corresponds to Pi{t) in the two-particle case and 
to Pi{t) in the single-particle case. 

We begin with the single-particle case for which, ac- 
cording to Eq.Q, the initial state (i.e. an eigenstate 
\Xi) of the isolated left well) can be expanded as 



■) = i^(o)) = E^"i^»'^)' 



(D2) 



where {\En')} are the eigenstates of the total system, 
i.e., of the single-particle Hamiltonian Hgp (pi). 



Next, we remark that our setup fulfills the so-called 
weak coupling condition between the left well (system) 
and the right well (environment): That is, for the low 
energies considered here, the projection of every eigen- 
state \En^^) onto the left well is either almost zero or 
approximately proportional to an eigenstate |x^) of the 
isolated left well. The latter implies that the level spac- 
ing A„ = e^''^^ — tn^'' of consecutive eigenstates of the 
isolated left well is much larger than the width of the 
corresponding Lorentzian distribution of the coefficients 
|c„p. Otherwise, the Lorentzian distributions of adja- 
cent low-lying levels would significantly overlap. In that 
case, the projection of an eigenstate \En^^) (with energy 
in the latter overlapping interval) onto the left well could 
yield superpositions of |Xj). 

We formulate the weak coupling condition for those 
eigenstates \En^'^) for which c„ ^ as 



E^^P\x)^bM{^)) ,forxe [0,£] 



(D3) 



hn being the proportionality factor. From |D3| and [D2| 
follows that hn approximately equals the expansion coef- 
ficient c„ 

Cn-hn- (D4) 

We highlight two major implications of the latter rela- 
tion: On the one hand, we find for each eigenstate \En^'^) 
of the total system with c„ ^ 0: 



Jo 



^(D5) 



If (as in our case) the coefficients |c„p are Lorentzian dis- 
tributed, so is the probability Pi{\En^^)) (for each non- 
vanishing component), see Footnote [7j 

On the other hand, we can show that Pi{t) and Psuiv{t) 
approximately coincide: 



pdt) = / dxm)\' 

Jo 

= Y^CnCn. rdxi?(-)*(x)i?L-)(x)e^(^^'-^^^'")* 
n,m 

fD5l 



E 



^ ^ |Cn| 6 I — -fsurv(0 i 



(D6) 



where we used that the {c„} and {&„} are real numbers, 
due to the real-valued Hgp in configuration space (for 
finite configuration space basis a real symmetric Hamil- 
tonian matrix) . The last line of ( D6 1 directly implies that 
(see, e.g. [T^) 

Psurv(t) = \TT[\cn\^6iE<^^P^ - 4^P))]|2 , (D7) 
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i.e., the survival probability is the (absolute value 
squared of the) Fourier transform (J-"T) of the so called 

local density of states pioc{E) = \cn\'^5{E^''P^ - E^f^^). A 
Lorentzian distribution of the coefhcients |c„p (see inset 
of Fig. [2| thus implies an exponential decay of Psurv (i) 
and — in the weak coupling approximation — an expo- 
nential decay of Pt{t). 



For two bosons, the identical calculation is performed 



by replacing tn' — 

V'fci ifc, 1 jx), and Pe{t) 



(sp) 



En, Xi{x) 



Pi{t). That is, the relations 
(|D5|) and (|D6|) directly carry over to the two-particle case, 



for which the survival probability is (almost) identical to 
the probability Pi {t) to find both bosons in the left well. 



